## ----setup, include=FALSE----------------------------------------------------------------------------------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
options(scipen = 999)


## ----eval = F----------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# if(!require("pacman", quietly = T)) install.packages("pacman")
# pacman::p_load(
#   # Data Manipulation & Cleaning
#   tidyverse, janitor, fs,
#
#   # Data Visualization & Plotting
#   ggthemes, ggridges, ggrepel, viridis, patchwork,
#
#   # Statistical Modeling & Analysis
#   brms, BradleyTerry2, icr,
#
#   # Utilities & Performance
#   remotes, tictoc
# )
#
# # Stan for bayesian models
# install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
# cmdstanr::install_cmdstan(overwrite = T)
#
# # Bayesian estimation for pairwise comparison
# remotes::install_github('davidissamattos/bpcs')


## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
pacman::p_load(
  # Data Manipulation & Cleaning
  tidyverse, janitor, fs,

  # Data Visualization & Plotting
  ggthemes, ggridges, ggrepel, viridis, patchwork,

  # Statistical Modeling & Analysis
  brms, BradleyTerry2, icr, bpcs,

  # Utilities & Performance
  remotes, tictoc
)

## Colors for visualization
ger_color  <- c(
  "AfD" = "#00ADEF",
  "CDU/CSU" = "grey60",
  "FDP" = "#ebcc34",# "#ffed00",
  "Greens" = "#46962b",
  "SPD" = "#E2001A",
  "Left" = "#8B1A1A"
)
# knitr::purl("README.Rmd", output = "main_script.R")

if(!dir.exists("data")) dir_create("data")
if(!dir.exists("model")) dir_create("model")
if(!dir.exists("plot")) dir_create("plot")



## ----eval = T----------------------------------------------------------------------------------------------------------------------------------------------------------------------------
model_dt <- read_rds("input/0_respondents_rating.rds") %>%
  glimpse


## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
mp_dt <- read_rds("input/mp_dt.rds") %>%
  glimpse()


## ----eval = T, include = T---------------------------------------------------------------------------------------------------------------------------------------------------------------
set.seed(888)
tic("Model training")
model <- bpc(
  data = model_dt,
  player0 = 'mp_1',
  player1 = 'mp_2',
  result_column = 'outcome',
  model_type = "davidson-U",
  cluster = "respondent",
  solve_ties = "none",
  priors = list(prior_lambda_std = 3.0,
                prior_S_std = 3.0,
                prior_U1_std = 3.0),
  show_chain_messages = T,
  seed = 888
)
toc()

write_rds(model, "model/0_main_model.rds")

tic("Parameter extraction")
param <- get_parameters_df(model, n_eff = T, Rhat = T)
toc()
write_rds(param, "data/0_main_params.rds")

main <- param %>%
  mutate(param_type = str_extract(Parameter, ".*?(?=\\[|$)"),
         page_id = str_extract(Parameter, "(?<=\\[A)\\d+"),
         respondent = str_extract(Parameter, "(?<=\\[A\\d{1,4}\\,)\\w+(?=\\])")) %>%
  glimpse()

lambdas <- main %>%
  filter(param_type == "lambda") %>%
  transmute(pageid = as.numeric(page_id),
            mean = -Mean,
            median = -Median,
            low = -HPD_lower,
            high = -HPD_higher) %>%
  left_join(mp_dt) %>%
  mutate(party = fct_reorder(party, mean)) %>%
  glimpse()

write_rds(lambdas, "data/0_mp_params.rds")

respondent_bias <- main %>%
  filter(param_type == "U1") %>%
  transmute(respondent,
            pageid = as.numeric(page_id),
            mean = -Mean,
            median = -Median,
            low = -HPD_lower,
            high = -HPD_higher) %>%
  left_join(mp_dt) %>%
  mutate(party = fct_reorder(party, mean)) %>%
  glimpse

write_rds(respondent_bias, "data/0_respondent_bias.rds")


## ----eval = T----------------------------------------------------------------------------------------------------------------------------------------------------------------------------
tic("Model BT")
model_bt <- bpc(
  data = filter(model_dt, outcome != 2),
  player0 = 'mp_1',
  player1 = 'mp_2',
  result_column = 'outcome',
  model_type = "bt-U",
  cluster = "respondent",
  solve_ties = "none",
  priors = list(prior_lambda_std = 3.0,
                prior_S_std = 3.0,
                prior_U1_std = 3.0),
  show_chain_messages = T,
  seed = 888
)
toc()
write_rds(model_bt, "model/1_model_bt.rds")

tic("Params BT")
param_bt <- get_parameters_df(model_bt, n_eff = T, Rhat = T)
toc()
write_rds(param_bt, "data/1_param_bt.rds")


## ----eval = T----------------------------------------------------------------------------------------------------------------------------------------------------------------------------
bt_data <- model_dt %>%
  select(mp_1, mp_2, outcome) %>%
  filter(outcome != 2) %>%
  mutate(outcome = ifelse(outcome == 0, "win", "loss")) %>%
  count(mp_1, mp_2, outcome)  %>%
  pivot_wider(id_cols = c(mp_1, mp_2), names_from = outcome, values_from = n)  %>%
  mutate(loss = ifelse(is.na(loss), 0, loss),
         win = ifelse(is.na(win), 0, win)) %>%
  arrange(mp_1) %>%
  mutate(mp_1 = as.factor(mp_1)) %>%
  arrange(mp_2) %>%
  mutate(mp_2 = as.factor(mp_2)) %>%
  glimpse

tic("Frequentist Bradley-Terry")
bt_model_obj <- BTm(cbind(win, loss), mp_1, mp_2, data = bt_data)
toc()
write_rds(bt_model_obj, "model/2_frequentist_bradley_terry.rds")

bt_preds <- BTabilities(bt_model_obj) %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  clean_names() %>%
  transmute(pageid = as.numeric(str_extract(rowname, "\\d+")),
            mean_bt = -ability,
            low_bt =  mean_bt - 1.96*s_e,
            high_bt = mean_bt + 1.96*s_e) %>%
  glimpse

write_rds(bt_preds, "data/2_frequentist_preds.rds")
